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Abstract : The paper describes a formulation of discrete scalar wave propagation in an 
inhomogeneous medium by the use of elementary processes obeying a discrete Huygens' prin- 
ciple and satisfying fundamental symmetries such as time-reversal, reciprocity and isotropy. 
Its novelty is the systematic derivation of a unified equation which, properly tuned by a 
single parameter, leads to either the Klein-Gordon equation or the Schrodinger equation. 
The generality of this method enables one to consider its extension to other types of discrete 
wave equations on any kind of discrete lattice. 

Resume : Cet article formule un modele discret de propagation d'ondes scalaires dans un 
milieu heterogene a l'aide de processus elementaires qui obeissent a un principe de Huygens 
discret et satisfont a certaines symetries fondamentales telles que le renversement temporel, la 
reciprocite et l'isotropie. Son originalite consiste en la derivation systematique d'une equation 
unifiee qui selon la valeur d'un seul parametre conduit soit a 1'equation de Klein-Gordon, soit 
a 1'equation de Schrodinger. La methode est suffisamment generale pour pouvoir envisager 
son extension a d'autres types d'equations d'onde sur toute sorte de reseau discret. 



I. INTRODUCTION 



The formulation of wave propagation by the use of Huygens' principle has stimulated extensive work in the 
past. Such an approach was first pioneered by P. B. Johns and coworkers who introduced the Transmission 
Line Matrix Modeling method (TLM) to solve the Maxwell equations in large electromagnetic structures 
j| which was later extended to the description of diffusive processes . Linear wave propagation is mediated 
by short voltage impulses which propagate along the bonds and are scattered on each node of a Cartesian 
mesh of transmission lines. The main idea behind such a formulation is a discrete equivalent of Huygens' 
principle, namely a principle of action-by-proximity obeyed by the pulses when they propagate at each time 
step from one node to the neighbor nodes and the emission of secondary wavelets at each node as described by 
the scattering process which radiates the incident energy in all directions. Besides the appealing viewpoint 
of Huygens' principle, the advantages of the method as a powerful tool for the numerical computation of 
wave propagation in complex and large structures have already been emphasized in the TLM method || . A 
similar approach in which the pulses were just scalar quantities propagating along the bonds of a Cartesian 
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lattice, was recently advocated for studying time-dependent wave propagation in large inhomogeneous media 
1-0- 

Here, we focus on the nature of the wave equations which result from such a formulation, especially on the 
possibility of retrieving standard wave equations. This question has already been raised in previous works 
where the scalar wave equation, the Klein-Gordon equation and Maxwell equations have been exhibited. In 
the TLM method, the equivalence between the voltage and current impulses with the electric and magnetic 
fields of Maxwell equations was initially demonstrated on a two-dimensional mesh |lj . The method was then 
extended to three dimensions and different features have been incorporated by many authors in order to 
describe losses, inhomogeneous media and boundary conditions [||. Another recent formulation || lead their 
authors to the classical wave equation, the Klein-Gordon equation and to the Schrodinger equation. However, 
due to restricting assumptions, these last results were limited to scalar waves in homogeneous media, and 
the Schrodinger equation was ill-defined at the continuum limit. 

This paper describes what is essentially a formulation from "first principles" of such an approach. The 
aim is to generalize the previous results to various kinds of waves in inhomogeneous media by considering 
the most basic properties obeyed by the scattering nodes. Such a formulation not only leads to the wave 
equations already mentioned but also results in a time-dependent Schrodinger equation which is well defined. 
Another advantage is the possible generalization of the method from the scalar waves considered in this paper 
to spinor || or vector waves. 

At each time step, the time-dependent discretized wave is defined as a linear combination of incident 
currents at each node of the Cartesian lattice. The nature of the wave depends on the number of distinct 
propagating currents on each bond of a given node. For example, one or several kinds of currents impinging 
on one node by the same bond will describe a scalar, a spinor or a vector field. Since they need only one 
type of current, scalar waves correspond to the simplest situation. For the sake of illustration, they will be 
the focus of this paper. In addition to the propagating currents, a current is attached to each node in order 
to describe an inhomogeneous medium. This current does not propagate but participates in the scattering 
process which involves the propagating currents. We shall demonstrate that a few elementary properties 
and symmetries of the scatterers like time reversal symmetry, reciprocity and isotropy lead to the classical 
wave, Klein-Gordon or Schrodinger equations. At that stage, we find that these discrete equations are not 
yet completely general because the mass, the velocity or the potential remain correlated. We show that it is 
possible to remove this limitation by attaching a second current to each node. 

The paper is organized as follows. In section [n| we introduce the currents which propagate on the mesh 
of a Cartesian lattice. Scatterers are described by scattering matrices located on each node of the lattice. 
The field is then defined as a function of incident currents. In order for the field to obey a discretized wave 



equation, the currents must be eliminated from its time evolution. We show in section III how this condition 
strongly compels the form of the scattering matrices. Time-reversal symmetry, reciprocity and isotropy of 
the scatterers are in trod uced in sections [fv|, |v| and VI respectively. The resulting discrete wave equatio ns are 



discussed in section VII 



This leads us to attach a second current to each node as described in section VIII 



We conclude in section IX by discussing possible extensions of this work. 



II. CURRENTS AND FIELDS 



A. Definition of the currents 



A current is a real or complex number, which propagates along a bond of a (i-dimensional Cartesian 
lattice. For scalar waves, it is sufficient to consider that each lattice bond carries two currents propagating 
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in opposite directions (Figure |l|).^| At any time t, the system is completely defined by the values of all 
currents in the lattice. We assume that time and space variables are discrete. In one time step r, a current 
propagates between two nodes of a given bond. All currents are synchronized. In other words, the currents 
are simultaneously outgoing from the nodes at some time t and become incident currents on neighbor nodes 
at time t + r. 



I 1 . 



MM 

FIG. 1. Currents propagating along the bonds of a two-dimensional Euclidean lattice. 



Each node of the lattice is a scatterer. The scatterers are identical in a periodic system or are distinct, and 
randomly chosen in a disordered quenched system. Each scatterer is described by a 2d x 2d matrix S which 
transforms the 2d incident currents Ei, at any time t, in 2d outgoing currents Sk at the same time (Figure 
||) according to 



2,1 



Sk = J~] su Ei k=l,...,2d, 



i=i 



where the Ski are the matrix elements of the matrix S. 
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FIG. 2. Scattering process. 



1 The results described in this paper are valid for any d-dimensional space. However, for convenience, all figures will 
represent a two-dimensional Cartesian lattice. 



3 



It turns out that the model is not sufficient to describe wave propagation in an inhomogeneous medium 
like a material in which the index of refraction is a function of the position. As shown later, this problem is 
solved by attaching to each node an additional current that will be called throughout this paper the node 
current. For notational convenience in sums over the currents, we shall write the node currents -f?2d+i and 
S^d+i- However, we emphasize that the subscript 2c? + 1 is not ascribed to a spatial dimension in contrast 
with the other indices ranging from 1 to 2d. The node current participates in the same scattering process as 
the other currents but does not propagate to neighbor nodes. Instead, the outgoing current S^d+i becomes 
the incident current E^.d+\ on the same node at the next time step. In other words, the propagation process 
for the node current reads E2d+i(t + t) = S^+iW- We can imagine the node current as propagating in one 
time step along a loop attached to the node. This loop is equivalent to the "permittivity stub" introduced 
in the TLM method to describe inhomogeneous electromagnetic structures [|10| . 

In conclusion, the general scheme of the time evolution of the currents is suitably represented in Fig. |^. 
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FIG. 3. Time evolution of the currents. 
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We write 

2(2+1 

S k =^s kl E u k = l,...,2d+l, (1) 

where the Ski are the matrix elements of the (2d + 1) x (2d + 1) scattering matrix S. This matrix is a function 
of the position in an inhomogeneous medium. 



B. Definition of the field 



The previous definition of the currents includes the three essential features of the model : the principle of 
action-by-proximity, which is described by the current propagation between neighbor nodes, the scattering 
process in which the nodes act like secondary sources of "spherical" wavelets and the linearity of this pro- 
cess. These features can be viewed as a discretized realization of Huygens' principle. However, additional 
ingredients and choices are needed to achieve the description of the model. For example, the form of the field 
propagation equation will strongly depend on the choice of the scattering matrices. Among all the questions 
to be answered in our approach, the first one is how to define the field. 

The field is supposed to obey a discretized linear wave equation which depends on the physical problem 
we are interested in : the scalar wave equation for an acoustic wave, the Klein-Gordon equation for a zero- 
spin massive particle, etc... At time t, the field is defined as a real or complex quantity -0(ii,Z2, ■ ■ • ,id,t) 
on each node of the lattice, where ii, ■ ■ ■ , id, are the discrete indices which locate the node position 
Xk = ikl, k = 1, . . . , d, as a function of the mesh parameter /. For briefness, we note r the set of indices 
ii, ii, ■ ■ ■ , id, of an arbitrary node. There is no reason for the field to be identified with one of the currents 
previously defined. Actually, we just postulate that the value of the field ip( r , t) 1S a function of the incident 
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and outgoing currents at node r and at time t. Therefore, by taking account of linearity, the most general 
definition reads 

2d+l 

i/>(r,t) = J2 [X' k (r)E k (r,t) + \Z(r)S k (r,t)], 

k=l 

where X' k and A' fc ' are complex coefficients to be determined. 

The field is a linear superposition of the currents, which vanishes when all currents vanish. By noticing 
that, due to the scattering process (Q), the outgoing currents Sk(r, t) can be expressed as linear combinations 
of the incident currents Ek(r,t) at the same time t, it is sufficient to define the field solely as a function of 
the incident currents 

2d+l 

i>(r,t)= J2 A fc (r)£ fc (r,f). (2) 
fc=i 

It is necessary to point out that the are functions of the node position. This property is required to 
describe an inhomogeneous system. 



III. CLOSURE OF THE WAVE EQUATION 

The general form of the discretized wave equations we are looking for can be written : 
4>(r, t + r) = f ty(r', t>), j,(r", t")^{r"', if"), . . .) , 



(3) 



where the field at node r and time t + r is a function of the field on the same node and/or on nodes at previous 
times t', t" , etc. . . Such an equation is straightforwardly derived from the discretization of the corresponding 
continuous equation. Equation (J3J) is a closed equation, which means that / is a function where no currents 
appear explicitely (so that (||) involves fields exclusively), and that the number of terms in the r.h.s. of (||) 
must be finite. Usually, and at the lowest order of the discretization procedure, the involved fields are fields 
on the same node r = ■ ■ ■ , id} and on the neighbor nodes i\ ± 1, i% ± 1, . . . , i<z± 1, at times t and t — t. 

In the remainder of this section we derive, by use of our model, a closed equation of the form of equation (^). 
Also, we show that imposing the closure conditions strongly constrains the form of the scattering matrix. 

For this purpose, we start to derive an equation of the type of Eq. (|[) from the evolution rules summarized 
in Fig. U of the model. First of all, let us introduce some notations. In a d-dimensional Euclidean space, node 
r is linked to the neighbor nodes by 2d bonds. We note the set of indices of the neighbor node which is 
linked to node r by the bond fc, along which the currents Ek(r) or Sk(r) propagate (Figure ^). Moreover, 
bond k of node r is named bond k for node rfc, which simply means, for example, that bond 2 of the central 
node in Figure |] corresponds to bond 1 of the right side node (k = 2,k = 1). This notation immediatly 
leads to Er(rk,t + r) = 5fc(r, t) and Ek(r, t + r) = Sj{rk,t). Note that for the node current (k — 2d + 1), 
we have r*2d+i = r and k = k. 



bond k 



bond k 



bond k 



propagation 



node r S k {r y t) S^r^t) node r 



bondl k 



node ri 



noder Ek {r, t + r) £f( r k , t + r 



time t 



time t + r 



5 



FIG. 4. Notation for the bond linking the central node r to one of its neighbor nodes r k and the propagating 
currents defined on this bond. 



The first step in attempting to derive an equation like ([|), is to write down the definition of the field 
ip(r,t + t) on a node r at time t + r 



2d+l 

Xk(r)E k {r,t + T), 


fe=i 




2d+l 




E Mr)Sj(r k ,t), 


fe=i 




2d+l 


"2d+l 


fc=i 


E s kA r k)Ei{r kl t) 
. 1=1 



(4) 

where the propagation and the scattering rules have been used. 

Thus, the field ip(r, t + r) at time t + r is expressed in terms of the incident currents at time t on the same 
node r and on its neighbor r k . One realizes immediately that using the same rules to express the incident 
currents at time f as a function of the incident currents at time t — t and iterating the procedure, will give 
incident currents at times t — r, t — 2t, . . ., on nodes as distant from node r as desired. Such a process is not 
bounded and seems to prevent us from writing a closed equation for the field ip. This statement is true except 
for some particular scattering matrices of the type we describe below. For this purpose, we shall rewrite the 
scattering equation ([!]) in order to exhibit the field ip : 

2d+l 2d+l 

k=i Am i=i Am 

where m is any of the 2d + I current indices. Suppose that the following condition is satisfied : 

= constant = p k , VI, (6) 

A/ 

then, equation (|J) becomes : 

S k = p k ip, 

and substitution in equation (Q) leads directly to : 

2d+l 

ip(r,t + T) = J2 A fc (r)^(r fe )i/;(r fe ,<), (7) 
fc=i 

which is a closed equation of the type we are looking for. It turns out that condition (||) is too restrictive 
and can be replaced by 

= constant = p k , VI ^ k, (8) 

Then, equation (||) becomes 

Sk = Pki> - PkE k , (9) 

where fik = PkXk - Skk- 

According to (||), the outgoing current on one bond is not only a function of the field i/j but also depends 
on the incident current on the same bond as sketched Figure ^| 
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FIG. 5. The scattering process depicted according to equation 



It seems at first sight that we went backward in writing condition (g) instead of (^) , since a current appears 
in equation (j^) in addition to the field ■0. Actually, substitution of equation (^) in equation (Q) gives 



2d+l 



2d+l 



(10) 



fc=i 



fc=i 



which is represented schematically on Figure 



2d+l 

ip(r, t + r) = x k{r)Pk( r k)^{rk, t) 
k=i 




k{r)ftj:{rk)Ej(r ky t} 
r r k 



FIG. 6. Sketch of equation (0). 

It is obvious from Figure [6] that the 2d+ 1 incident currents Ej(rk, t) at time t (bold arrows) are deduced, 
according to the propagation rules, from the currents Sk(r,t — r) leaving node r at time t— r, as depicted 
in Figure 




FIG. 7. The outgoing currents leaving node r at time t — r, which give after propagation the input currents at time 
t in equation (J10|) . 
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Thus, the additional currents due to the looser condition (|S|) eventually involve onlythe initial node r 
instead of propagating new terms to remote nodes. Using again equation (|9|), equation ( |10| ) becomes 



2d+l 



ip(r,t + r) = x k{r)pj;(r k )ip(r k ,t) 



k=l 



'2d+l 



£ x k(r)p k (r k )p k (r) 



k=l 



2d+l 



ij){r,t-r)+ Y ^k(r)p k (r k )p k (r)E k (r,t - r), 



(11) 



k=l 



The last term is a linear superposition of 2d + 1 incident currents upon node r at time t — r. In order for 
equation (Jill) to be a closed field equation without current terms, it is sufficient for this superposition of 



currents to be proportional to tp(r, t — r) = Y^k=i ^k(f)E k (r, t — t). Therefore, we must set 

K {r k )^ k (r) = fj, 2 (r), Vk = 1, . . . , 2d + 1, 



(12) 



where we denote p 2 (r) a constant which is a priori a function of r. In fact, it is simple to prove that p 2 (r) 
does not depend on r. To prove this we begin by recognizing that equation ( [l2| ) can be written for one of 
the neighbor node r k 



Hj{r k i)ni{r k ) = pL 2 (r k ) 



V/ 



,2d+ 1. 



Now r k plays the role of the central node r of equation (|lj) and r k i denotes its neighbor nodes. Among the 
neighbor nodes r k i, the initial node r (7 = k) is of particular interest and the last equation reads 

Vk{r)^ k {r k ) = li 2 {r k ). 

Comparison with equation ( [l2|) leads to the required identity 

M 2 (r fe ) = M 2 (r). 

Since the node r is arbitrary, equation ( |l2| ) can be recast in a simpler form valid for any node 

Hk{rk)Hk{r) = [i 2 , Vfc, 



(13) 



where now the constant /j, 2 is independent of the node's position. 

One notices that for k = 2d + 1, = r, k = k, and ^ k {r k ) — ^ k (r). Therefore 

V2d+i(r) = e 2 d+i(r)v, 

where e 2 d+i(r) = ±1. 

Finally, utilizing equation (|l3|), equation ([ll]) reads 



(14) 



2d+l 



i/j(r, t + r)= Y ^k(r)pj:{r k )iP(r k , t) + /i 2 



fc=i 



2d+l 

£ 

k=l 



Xk(r)p k (r) 
Mr) 



ip{r, t — r). 



(15) 



Equation ( [l5| ) is a closed field equation in which currents do not appear. As the field ip is exhibited at 
times t + r, t and t — r, this new equation is an improvement over equation (R) in the sense that a discrete 
second time derivative ijj(t + r) — 2ip(t) + ip(t — r) can be obtained as required in, for example, the discrete 
time-dependent wave equation. Making use of condition (ra) , the general matrix element reads 



Ski — Pk^l ~ Mfc^fci, 



(16) 



where = p k \ k — s kk and 5 k i is the Kronecker symbol. Instead of (2d+ l) 2 elements, the scattering matrix 
depends now on the 3(2d+ 1) complex coefficients p k , X k and fi k . These coefficients vary independently from 
one node to another in an inhomogeneous system except for the coefficients p, k which obey condition ( |l3| ) 
relating neighbor nodes. Furthermore, pb 2 d+i = ±Mi according to equation (H). 
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In summary, equation (|15|) is the general equation which describes the time evolution of the field ip. This 
closed equation has been derived from the evolution rules of the currents, the definition of ip (Q) and from 
conditions (||) and ([l3]). Let us discuss the status of these conditions. Condition (||) cannot be avoided if we 
look for a closed equation in ip, since a looser condition would make terms appear on all nodes of the system. 
Meanwhile, condition (|L|), which allows one to replace the currents in the last term of equation (^lj) by the 
field ip(r,t — r), is sufficient but not necessary. Indeed, suppose we consider that these incident currents 
defined on node r at time t — r are the result of propagation and scattering of currents at previous times 
t — 2t and t — 3r on the neighbor nodes and on node r. Reproducing the same calculation steps as those 
which lead to equation (|ll]) , the last term may be replaced on the one hand by new terms involving the field 
ip on node r at time t — 3r and on its neighbors at time t — 2r respectively, and on the other hand by a new 
linear superposition of incident currents on node r at time t — 3r. Therefore, it is possible to recover a new 
closure condition (analogous to equation (|l3|)) by taking this new linear superposition of incident currents to 
be proportional to ip(r,t — 3r). Additionaly, the same procedure can be iterated and similar superpositions 
appear at times t — 5r, t — It, etc. . . So that the resulting field equation can be closed at any step of the 
iteration procedure. The net result is to open the possibility of describing high order time derivatives of the 
field ip. We shall not consider these possibilities further and restrict ourselves with equation ( Jl5| ) which is 
sufficient to obtain the second order time derivative of ip. 

The description of the model is far from being achieved. To make further progress, additional properties 
of the scatterers are needed. The most basic properties to consider are suitable symmetries of the scattering 
process which will restrict again the number of independent parameters and make precise the physical content 
of equation ( |l5| ) . We shall successively introduce some of the most natural symmetries we can think of, namely 
time-reversal symmetry, reciprocity and isotropy. It is worthwile to point out that such symmetries, like the 
results obtained so far, will only involve the discrete geometry of the Cartesian lattice and assume neither 
any topological space nor any Euclidean or Riemannian structure. Nevertheless, we shall need the definition 
of underlying manifolds provided with such structures when taking later the continuous limit of the discrete 
wave equation obeyed by the field ip. 

IV. TIME-REVERSAL SYMMETRY 

In this section, we introduce the time-reversal symmetry by formulating its natural action on both currents 
and fields. 



A. Currents 

Time-reversal is naturally completed by reversing the direction of the current arrows in such a way that 
the currents propagate and are scattered backward in time. They will obey time-reversal symmetry only if 
they describe the past states of the system in reverse order. It is then obvious that the scattering process 
must be reversible. Therefore, we state that the forward scattering process depicted in Figure || (a) implies 
the existence of the reverse process depicted in Figure || (b). Since the currents are, in general, complex 
quantities, complex conjugate currents are involved in the time-reversed process. If the currents are real, this 
will be equivalent to reversing the original currents. 
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FIG. 8. Processes involved in the time-reversal symmetry, (a) forward process, (b) time-reversed process. 



Using matrix notation, time-reversal invariance reads 

(S k ) = S(E k )^(E* k )=S(S* k ), 

where (S k ) and (E k ) are the column vectors of outgoing and incident currents respectively. 
Then, taking the complex conjugate (Ef.) — S* (S k ), one obtains 



SS* = S*S = I, 



2d+l, 



where 22<i+i is the (2d+ 1) x (2d + 1) unit matrix. 

This condition and the general form ([l6|) for the scattering matrix elements leads to the following conditions 



Pk 



= constant = e* 71 , 



constant 



,2d+l, 



(17) 



and 



2d+l 



>*7i 



"272 



(18) 



fc=i 



Again, let us point out that all the quantities appearing in these relations are functions of the node r. This 
statement holds in particular for the constants C\ = e 171 and C2 = e 172 . 
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B. Time- reversal symmetry of the field 



Equation (plf ) which describes the time evolution of the field can be viewed as some general relation which 
expresses the field ip at time t + r as a function of ip at times t and t — r 



iP(r,t + t) = f [^j(r k ,t),iP(r,t - t)} 



(19) 



In general, ip is complex valued function so that equation ( |19| ) obeys time-reversal symmetry if the reverse- 
time equation is verified by the complex conjugate field ip*, i.e., 



if (r, t - r) = / [V* (r fe , i), V* (r , t + r)] 



(20) 



Therefore, our task is to check that the time-reversal symmetry translated on currents, which was considered 
in the previous section, enables us to obtain equation ( pp| ) from equation (|l9|). For this purpose, we first 
substitute the definition (|J) of the field ip in ( |l9|) 



2d+l 



53 X k (r)E k (r,t + r) = f 



fc=i 



'2d+l 



2d+l 



53 k(r k )Ei{r k ,t), J2 Xk(r)E k {r,t- T ) 



i=i 



fc=i 



(21) 



Then, Figure || tells us that reversing the direction of current propagation and substituting the complex 
conjugate currents S k for E k gives the current values of the reverse-time system. In other words, the time- 
reversal symmetry applied to the currents, appearing in equation (^l|), leads to 



2d+l 



53 X k (r)S*(r,t-T) = f 



k=l 



2d+l 



2d+l 



53 h(r k )St(r k ,t), J2 \k(r)S*(r,t + T) 



i=i 



fc=i 



(22) 



This equation is similar to equation (pp|). However, it involves the quantities X^fefii 1 X k S k instead of ip* = 
Efcfii 1 ^*k^k- Therefore, to demonstrate the time-reversal invariance of the field ip, it is necessary to show 
that these two quantities are identical. Actually, they just need to be proportional to each other since the 
function / is linear. Using the scattering formula (^|), we can write 



2d+l 



2d+l 



53 x k s* k = A,r - 53 \ k »* k E* k 



k=l 



k=l 



Since X k /J, k = C|A£ according to the second equation of ([n]), one obtains 



2d+l 



53 X k S* k = (Ai - C 2 *)V* = dp , 



k=l 



where we have used equation (|T^). However, we have previously pointed out that C± is a function of the 
node r. Therefore, substituting C\ip* for Y^ 2 k =i XkS^ in equation (|2^ ) leads to 

& (r)r (r,t-r) = f[C 1 (r k )ip* (r k , t), & (r)iP* (r , t + r)] . 

This equation will be equivalent to equation ( p0| ) if we impose the condition 

Cx{r k ) = d(r), Vfc. 

Thus, the constant C\ must be the same on the node r and on its neighbors r k , for any node r in the system. 
One must conclude that C\ is a constant over the whole lattice. Therefore, the first equation of ( |l7| ) becomes 



Vk(r)p* k (r) 
Pk{r) 



Vr,k. 



(23) 



11 



C. Form of the time-dependent equation 



We now proceed to show that the results obtained so far have already strong consequences on the form of 
the evolution equation of ip. First, the coefficients of ij ){r, t — r) and if)(r k ,t) in equation ( fl5| ) can be recast 
in the following form, with the help of identities (|l7|), (18) and (53) 



2d+l 



Afc(r-) / Qfc(r) 
MfcO") 



Ci<7 2 (r) 



and 



Afc(r)pfc(r fe ) = e fc (i-)|Afc(i-)||^(rfc)| 



CiC 2 (r) 



1/2 



where the modulus and phases of A and p have been separated in the above expression. The sign e k (r) = ±1 
depends on the choice made in taking the square root. Plugging these two results into ( p!s| ) leads to the 
time-reversal invariant equation of the field 



2p(r,t + T) + 



2d+l 



CiC 2 (r) 



V>(r,t-r) = x k{r)pj{r k )t)}{r k ,t), 



k=l 



CiC 2 (r) 



1/2 2 ci+l 

fc=i 



(24) 



where |/i 2 /(CiC 2 (r))| = 1 since |^| = |d| = \C 2 (r)\ = 1. 

It is now easy to convince ourselves that the form of ( pi[ ) is either d 2 Tp{r, t) / dt 2 = Lip(r, t) or idip(r, t) jdt = 
Lip(r, t), where L is a real operator. Let us notice that the left hand side of equation (24) contains the terms 
which arc needed to build the time derivatives of the field. There are only two possibilities : p 2 / [CiC 2 (r)] = 1 
or I [CiC 2 (r)] = —1. If I [CiC 2 (r)] = 1, the left hand side of equation (|24|) becomes ip(r, t+T)+ip(r, t—r) 
which allows us to construct a second order time derivative : [ip(r, t + t) + tp(r, t — r) — 2ip(r, t)] /t 2 . If 
A* 2 / [CiC 2 ( r )] = — 1, the difference tp(r,t + T) — i/j(r,t — T) is directly related to the first order time derivative : 
[ip(r, t + r) — ip(r, t — r)] /2r. The first or second order time derivatives of the field must be the same on the 
whole lattice, which is achieved by setting C 2 (r) = C 2 . 

Lastly, we conclude that according to the two values we assign to /Li 2 /(CiC 2 ), equation fl24|) becomes 



(CiC 2 ) 



2d+l 



ip(r,t + r) + ip(r,t— r) 



k=l 



(CiC 2 ) 



2d+l 



i[^(r,t + T)-V(r,t-r)] = ^2e k (r)\X k {r)\\p^{r k )\i;{r k ,t), 



k=l 



as was claimed above. In particular, it is worthwhile to point out that, due to the square root of /x 2 /(CiC 2 ) 
in the right hand side of (pjj), the factor i appears when the equation is of first order in time. 
In the remainder of the paper, we shall write 



(CiC 2 ) 



^ 2e -i( 7l + 73 ) = £e 



(25) 



with e e = ±1. 
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V. RECIPROCITY 



In this section, we go further by imposing a new symmetry, reciprocity, on the evolution rules of the 
currents. The definition is sketched in Figure ^ where E and S are respectively an input and an outgoing 
current defined on two arbitrary bonds of the same node. We assume that the direct process in Fig. KJ (a) 
implies that the reciprocal process in Fig. || (b) also exists. In the reciprocal process, the currents E and S 
have the same values as in Figure ^ (a) but the two arbitrary bonds have been exchanged. 




scattering 



s 



t 



s 



/ 



+ 



(a) 



(b) 



FIG. 9. Reciprocity, (a) direct process, (b) reciprocal process. 



Formulated in this way, this symmetry is the counterpart for the scattering process, of the reciprocity 
theorem which is well known for propagation in inhomogeneous media (see for instance JTT[-^| ) . This is our 
main motivation for introducing it. Let us discuss its consequences. 

Using matrix notation, the process described in Fig. ^ (a) can be written as Si — sn~Ef. where E^ is the 
input current and Si is any of the outgoing currents. We can also imagine a second process for which the 
only non-vanishing input current is incident on bond I while we are looking at the outgoing current on bond 
k, i.e. Sk = SkiEi. If both processes are chosen so that Ek = Ei, reciprocity tells us that Sk — Si, leading to 
sm = sik- Since the two bonds k and I are arbitrary, one concludes that the scattering matrix is symmetrical 

3 = S T . 

Additionaly, the scattering matrix is time-reversal invariant : S^ 1 — S* . Hence, the scattering matrix is 
unitary. As a direct consequence of this result, the sum of current intensities is conserved in any scattering 
process on the lattice 

2d+l 2d+l 

]T |^.| 2 = ]T \E k \ 2 . (26) 

k—l k=l 
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Since the current intensities are also conserved in the propagation stage, their total sum over the whole lattice 
keeps the same value at each time step. This property guarantees the stability of the current dynamics. To 
understand its importance, let us consider a system of linear size L. The linear evolution of the whole set of 
currents can be described by a (2d + l)L d x (2d + l)L d matrix which transforms the (2d + l)L d input currents 
at time t into the (2d + l)L d input currents at the next time step t + r. Such a linear system is characterized 
by its eigenvectors and the corresponding eigenvalues. Any arbitrary initial state of the currents can be 
expanded on the basis of these eigenvectors. Let us choose one arbitrary eigenvector, characterized by its 
eigenvalue £> a s the initial state. The current intensities will diverge to infinity if | CI > 1 or decay to zero 
if | CI < 1- The above result, equation (^6|), prevents us from such a scenario since the total sum of current 
intensities is constant regardless of the initial state. We conclude that all the eigenvalues are bound to be of 
modulus 1 and that any initial linear superposition of eigenvectors is preserved by the time evolution of the 
currents. This result is pleasing because it is equivalent to what is expected from the acceptable solutions 
of the usual linear wave equations such as the classical wave or Schrodinger equations. In order to break 
the stability of the current dynamics, losses or interactions of the system with the "outside world" must be 
introduced in some way. 

Note that instead of introducing reciprocity, we could have required first the conservation law (|26|), and 
then combining this law with time-reversal symmetry would have led to reciprocity. 
From Ski = sik and (^), we find 

Pkh = PiAfe, Vfc, I, 



pk 

A A' 



i?e i7 , 



where R and 7 are independent of k. 

Combining this result with the two first equations of ([Ft]) leads directly to 2j = 72 — 71 . From ( |l8| ) and 
the definition (f25f) of e e , one obtains 



-H2 1 p 1 ~1 



R = 



Pk = 



f l + e e /p 2 ) p k X* k 
Z^fc=l |Afe| 



(27) 



The matrix element becomes 



Ski — PfcA/ — PkSki — Pk 



(l + e e //i 2 ) A£A, 
l^k=l lAfcl 



Remembering that, according to the second equation of (|l7]), pk = e472 Afc/Afc, the scattering matrix is 
parametrized by two constant phases 71 and 72 and by 2d + 1 complex coefficients Xk which are functions of 
the location r. 



VI. ISOTROPY 



Isotropy is not required in order for the model to be solvable. However, we impose isotropy to simplify 
the model since anisotropic scatterers could be considered as well. Let us consider a process where only one 
of the propagating input currents Ei is different from zero (Figure [l^ (a)). This process is by definition 
isotropic if all outgoing currents have the same value except for two outgoing currents : Si on the input bond 
and the node current Szd+x 

Sj = S k , Vj, k I j,k^i, j, k^2d+l. 
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In other words, when viewed from one particular bond, all other bonds are equivalent except the node bond 
because of its special status. The scatterer will be isotropic if the scattering process is isotropic for any chosen 
input bond. 

Another process to consider is the one where the node current is different from zero (Figure 0(b)). This 
process is isotropic if all outgoing currents have the same value except the node current S^d+i, i- e - 

Sj = S k , Vj,k I j,k^2d+l 

Actually, both kinds of processes lead to the same conclusions because of the form of the scattering elements 
discussed in the previous section. 



Pj = Pk, Aj = A fc , fij = /i fc , Vj, k I j, k ^ 2d + 1. 



(28) 



1 



(ij + l) 



> 







scattering 



t 
1 



(« + l,i) 





q scattering 



+ 



Si 



Si 



Si 



/ 

S2&+1 



Si 



t 

i 



Sk 
Sk 



Sk 



/ 



Sk 



t 



(a) 



(b) 



FIG. 10. The two scattering processes considered to define isotropy. 



The scattering matrix now depends only on two complex parameters Ai and \2d+i, and on the two phases 
71 and 72 (since ^ and pt can be expressed as functions of Xk according to the second equations of (njh and 

Conditions (EH) and (Q) provide an additional constraint for the value of p,\. First, suppose the value of 
pi is known at node r. Then, according to (|l3|), the value of /j,-^ on the nearest-neighbors of node r reads 



Pkirk) = 



Vfc ^ 2d+ 1. 



Pk(r) MiW 

Since the scatterers r& are also isotropic, /Ujr(rfc) = Mi( r fc) f° r k ^ 2d = 1, we obtain 



Vfc ^ 2d+l. 
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Thus, /Lti(rfc) has a constant value on the nearest neighbors r& of node r. Next, starting from nodes r*., we 
obtain for the second nearest-neighbors ru of node r, 

Hi{r u ) = — — - = /ii(r). 

Finally, proceeding in the same way from node to node, we find that the whole system is made of two 
intermingled Cartesian lattices characterized by their respective values pL\ and p 2 / p,\. It is important to note 
that this result is peculiar to a Cartesian lattice, which has been considered from the beginning. Actually, 
the model we have discussed so far does not depend on the type of the discrete lattice. Indeed, if we had 
considered, for instance, a triangular lattice, a quasiperiodic lattice or even a random lattice (where the node 
connectivity changes from node to node, see Appendix), the time dependent equation of the field or the form 
of the scattering matrices would not have changed. However, the conclusion leading to the two values of 
on a Cartesian lattice would need to be modified on an arbitrary lattice. Consider a triangular lattice in 
which the elementary cell contains three nodes connected by three bonds. Following an elementary closed 
loop of three bonds starting at node r, we would find after a complete cycle 

i \ V 2 
/W) = — 7-7, 



pi{r) = 



(29) 



where t\ = ±1. Moreover, e\ should have the same value at the three nodes and as a result over the whole 
system. As this conclusion is true for almost any arbitrary lattice, we are led to adopt ( p9| ) for the sake of 
generality. Remembering that p is defined by its square p 2 in equation (|l3|), it is possible to choose once and 
for all e% = 1, i.e. p\ = p 

Let us combine now these results with time reversal symmetry and reciprocity. For this purpose, we 
introduce the ratio v = A 2( j+i/Ai. From the second equation of (17) and (O), we write 



Mi = P2d+1 — = £2d+lM — 
v v 



Since p\ = p, we get 



— — £2(2+1 — ±1- 
v 

That is, the ratio v = A2d+i/Ai is either real or imaginary. 

On the other hand, reciprocity leads to p2d+i = P\v- Hence, /02d+iAi 
p\\\v 2 . Therefore, the scattering matrix elements become 



PiA 2 d+i = PiAi^ and p 2 d+iA 2 d+i 



sh = PiAi - pSki, 



\/k,ljt 2d+l, 



Ski = sik = pi\iv, 



Wk^2d+l,l = 2d+l, 



(30) 



where, using ( |27[) 



Ski = piXiv 2 ~ pe 2 d+i, k = l = 2d+l, 



PiAi 



2d+ vv* ' 



(31) 



Let us make a few remarks about these last results. First, the coefficients Xk which appear in the expression 
of the field ip(r,t) — Y^k=i Xk(r)Ek{r,t) are defined modulo a multiplicative constant. Obviously, the 
scattering matrix cannot depend on such an arbitrary constant, and indeed, one can see that it is only 



1G 



function of the ratio v = Xzd+x/Xx, according to ( |30| ) and Next, /i and e e are constants denned over 

the whole system. Hence, v is the sole quantity which depends upon the location r in the system (with the 
restriction v* jv = £2d+i = ±1)- Given the definition v = A2d+i/Ai, we see then that this is the node current 
which enables us to describe an inhomogeneous system as asserted before. A model without node currents, 
which corresponds to v = 0, would have restrained us to an homogeneous system. We can interpret the node 
current as trapping a fraction of the wave energy before releasing it to neighboring nodes. The parameter 
v{v) measures the local strength of this trapping effect. We shall see in the next section that the velocity of 
the wave is directly related to v(r). 

We can now write the final form of the wave equation ( fji]) as 



ip(r,t + t) + e e ip(r,t- r) = Ai(r) 



2d 



Pi( r k)^P(r k ,t) + v 2 (r) Pl {r)tj){r 1 t) 



Lfc=l 



Given this equation, we are led to introduce the new field 



4>(r,t) = px{r)i>{r,t) = pi(r)Xi(r) 
which obeys the following equation 



2d 



J2E k (r,t) + v(r)E 2d+ i(r,t) 



(r,t + T) + e e (/>{r,t-T) 



P + e e/p 

2d + v(r)v*(r) 



2d 



5^0(r fc> t)+i/ a (r)0(r,t) 



k=l 



Like the scattering matrix elements (30), this equation depends only on the constants fi and e e and on 
the local parameter v(r). We note that the coefficient Ai(r) comes into the expressions of the field, of the 
scattering matrix and of the wave equation as the product Ai(r)pi(r). Therefore, it is possible and convenient 
to rename the quantities \\{r)pi(r) and X2d+l(f)pi(r) as Ai(r) and A2d+i(^) respectively. With this new 
notation, we can reformulate the above equations as 



(r,t) = A 1 (r) 



2d 



Y J E k (r,t) + v{r)E 2d+1 {r,t) 



k=l 



Mr) 



H + e e /p 
2d + v(r)v*(r)' 



(32) 
(33) 



(r,t + r) + e e cf){r,t-T) = A x (r) 



2d 

^(0(r fe ,i)-0(r,t)) + (2d + ^(r)) ( /»(r,t) , (34) 

,k=l 

where we have shown explicitely the term X)fe=i ( ( / , ( T *fc! *) — <M r 7 1))> which corresponds to the discrete Lapla- 
cian on an Euclidean lattice. Similarly, the scattering matrix is given by 



5(r) = Ai(r) 



/ 1 



1 



1 v(r) \ 



1 v{r) 
v(r) v 2 (r) J 



(I 

o '•. 



\ 
: 



e 2 d+i(r) / 



(35) 



where £2d+i{ r ) = v*(r)/v(r) = ±1. 
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VII. DISCUSSION OF THE WAVE EQUATION 



The wave equation ( |34| ) is a function of two discrete parameters, e e and C2d+i(v). The constant e e is defined 
over the whole system. We have already seen that its value ±1 controls the order of the time dependent 
wave equation, which will later be discussed separately. The parameter e2d+i(r) has a different satus. Being 
defined by e2d+i{r) = v*(r)/v{r) with v{r) controlling the disorder of the system, its value ±1 depends 
on the node r. However, the type of disorder we would like to describe should depend continuously on the 
location in the system, such as a potential, a wave velocity or an index of refraction which varies smoothly 
with the position. Therefore, we shall assume that the value of e^d+i is fixed over the whole system whereas 
v{r) is a smooth function of r. If £2rf+i = 1 or —1, v(r) is real or imaginary, respectively. Both cases will also 
be considered separately. Hence, four possibilities will be discussed according to the values of e e and e2<j+i. 



A. Second order time equation (e e = 1) 



1 



Since in this case v(r) is real, v(r)v*{r) — v 2 (r). The wave equation becomes 

2 cos 9 



(r.t + T) 



(r,t- t) - 2(f>(r,t) 



2d- 



2d 

E 

k=l 



(<f>(r k ,t) - 0(r,t)) + 2(cos0 - l)0(r,i), 



where we have written explicitly the second order time derivative and fi = e zS . 

We have already pointed out that our results were not peculiar to a Cartesian lattice. A similar wave 
equation would have been obtained in any discrete lattice. However, in order to consider this equation as 
the discretized version of the corresponding continuous wave equation, we need now to limit ourselves to a 
Cartesian lattice so that we can introduce the lattice constant I and the current velocity cq = 1/t. Hence, 
one writes 

2d 



2 cost 



2d + v 2 (r) 
1 



§)Ew(^,t)-^,t)) 



fe=i 



I 2 



2(1 - cos( 



(r,t), 



(cf>(r, t + t) + cf>(r, t - t) - 20(r, t)) = 
which is the discretized version of the Klein-Gordon equation 

In the particular case where the Klein-Gordon equation describes a scalar particle of mass m, the coefficient 
a would be given by a 2 — m 2 c 4 /h 2 . Comparison between ( |36|) and ( |37| ) provides the velocity c of the wave 
and the coefficient a 



(36) 



(37) 



c 2 (r) 




(38) 



The special choice 9 = 0, i.e. /1 = 1, leads to the scalar wave equation c 2 {r)V 2 <t> — d 2 cf>/dt 2 = 0. The first 
equation of (|38| ) shows clearly that the velocity of the wave is directly connected to the parameter v(r). Its 
maximum value c max = co(cos9/d) 1 / 2 is reached when v(r) = 0. As stated before, the role of the node 
current is to slow down the wave. If d > 1, we also note that c max is always strictly smaller than the velocity 
Cq of the currents. 
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2. e2d+i = 



-1 



In this case v{r) is imaginary, thus v 2 (r) = — \v(r) \ . The wave equation becomes 

/ 2 cosO 9 \ ^ 1 , . , 

•^ Ep(^.t)-^.i)) 



v 2d+|i/(r)| 
1 



fe=i 



(0(r,* + r)+0(r,*-r)-20(r,i)) 



cos 6* 



2d+\u(r)f 



1 



(39) 



This equation is also a discretized version of ( p7| ) since the factor of <fi(r,t) in the right hand side of (|39|) is 
positive. The velocity of the wave is given by 



c 2 (r) 



2cos6> 



2d+\v{r)Y 



However, ( |39D is not as versatile as (36). First, in order to obtain the scalar wave equation, we must set the 
double condition (9 = and v{r) = 0. Hence, the wave velocity is pinned to its maximum value c = Co/"\/d. 
This describes a uniform system without node current. Next, it is easy to show that the same double condition 
holds if we take the continuous limit ( |39] ) in fixing the value of a 2 . Under these circumstances, we conclude 
that an imaginary node current is not compatible with equation (^7|). 

Hence, from the analysis of these two cases (e2d+i = ±1), m order to correctly describe the Klein-Gordon 
equation and the scalar wave equation in a inhomogeneous system by our current model, we must set e e = 1 
and e-2d+i = 1- The resulting equation is given by (|36]). This result is not completely satisfactory since the 
necessity of setting 62^+1 = 1 is not derived from the g eneral kind of reasoning which has been used in our 
approach up to this stage. We will see in section VIII that adding a second node current to each scatterer 
will remove this arbitrary condition. 



B. First order time equation (e e = — 1) 



1- £2d+l = 1 

With v(r) being real, the wave equation becomes 

- (0(r , t + r) - 4>{r, t - r)) = -- [^—^ ) £ [ p J - — <t>(r, t). 

k—1 

We recognize the discretized version of the Schrodinger equation 
with 

V{r) _ smO 

and 

h I 2 ( sin \ l 2 V(r) 
2^~~ \2d+\ V {r)\ 2 ) ~ ~ h{2d+\ V (r)\ 2 )' 

The expression for V(r) shows that we would obtain a Schrodinger equation in which the potential is bound 
to stay constant over the whole system. Therefore, we shall reject the choice e e = —62^+1 = — 1- 
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2. e 2d +i = -1 



Again, v 2 (r) = — \u(r)\ . We obtain 



sinfl { (f>(r k ,t) - <f>(r,t)\ sin(9 / 2d - | v{r) | 2 " 



^ n , fc , t ;-^, t; _ ™ ^~in^l 2 ^ (rjt)< (41) 



2d+|^(r)| 2 / ^ V i 2 / t~ \2d+\v{r)\ 

In a similar way, identification with the discrete Schrodinger equation (BO) gives 



V(r) _ sin 6 / 2d - 



2 



2d+|Kr)| 2 7 ' 



and 




l 2 V(r) 



h(2d-\v(r)\ 2 ) 



In contrast with the previous section, the potentiel V(r) is a function of the position r. The last equality 
shows that the mass m also depends on r. This is not a surprise in the context of quantum waves on a lattice 
where one usually introduces an effective mass which may be a function of position (see |l5[] for instance). 
However, the Schrodinger equation ( f4l| ) is not completely general. We note that the variations of m and V 
with r are correlated. In particular, the mass cannot be kept fixed at a constant value if the potential is a 
function of r. In the context of the localization of quantum waves in an inhomogeneous medium, equation 
(fll]) would correspond to the Anderson model with correlated "diagonal" and "off-diagonal" disorder fl5|| . 
Hence, one may wonder which conditions must be imposed in our approach in order to derive a Schrodinger 
equation where the two types of disorder can be described separately. The most natural idea which comes 
to mind, is that one needs an additionnal degree of freedom which would allow to separate the mass and the 
potential spatial variations. Introducing a s econd node current not only achieves this task but also removes 



the arbitrary condition discussed in section VII A 



VIII. INTRODUCING A SECOND NODE CURRENT 



Let us attach to each node additional currents E^d+i and 52(2+2 which behave like the node currents i?2d+i 
and 52ci+i- They participate in the same scattering process as the other currents but do not propagate to 
neighbor nodes. Instead, the outgoing current S2d+i becomes the incident current on the same node 

at the next time step. The scattering process and the definition of the field become 



2(2+2 



S k =J2 s kiE h k=l,...,2d+2, 

i=i 

2(2+2 

1&(r,t) = *k(r)E k (r). 



k=l 



It is easy to check that the derivation of a wave equation similar to (34) remains quite unchanged. Therefore, 
we will only point out the most noticeable differences and the final result. 
First, equation (|l4|) is supplemented by 

H2d+2(r) = e 2 d+2(r)At, 



where e 2 d+2(^) = ±1. 
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Second, instead of the single ratio v = Xad+i/^i, we introduce V2d+i = ^2d+i/^i and V2d+2 = ^2d+2/^i 
which obey the conditions u^d+i/^d+i = £2d+i = ±1 and v^ d+2 l v 2d+2 = ^2d+2 = ±1- In other words, we 
find again that V2d+i and V2d+2 are either real or imaginary. 

Finally, equations (B2), (B3T) and (134) become 



0(r,t) = Ai(r) 
Ai(r) 



' 2d 

E 

.fc=i 



E k (r,t) + v 2 d+i{r)E 2 d+i{r,t) + V2d+2{r)E 2 d+2{r , t) 
fi + e e /fi 



2d + \v2d+i(r)\ + \v2d+2(r)\ 



(r,t + T) + e e (f>{r,t~T) 



Mr) 



2i! 



v 2d+ii r ) + v 2d+2{ r )) <l>(r,t) 



.k=l 



(42) 
(43) 

(44) 



At this stage, the discussion of section VII can be reproduced. However, as shown below, we need only to 
discuss the two cases e e = ±1 separately. Indeed, we have been previously led to discuss the two choices 
v = A2d+i/Ai being real or imaginary. By adding a second node current, it happens now that we have gained 
additional possibilities since i^d+i = A2d+i/Ai and V2d+2 = A2d+2/Ai can be real or imaginary independently 
of each other. However, one can easily convince oneself that choosing ^2rf+i and V2d+2 being simultaneously 
real or imaginary is equivalent to having a single node current such that v 2 = v^d+i + v 2d+2- The only new 
possibility opened up by introducing two nodes currents corresponds to V2d+i real and ^2d+2 imaginary, or 
vice-versa. Since both node currents have been equivalent until now, it is sufficient to consider one of both 
cases, i.e v 2 d+i real (e 2 d+i(r) = 1) and V2d+2 imaginary {e 2 d+2{r) = -1). 



A. Second order time equation (e e 



By introducing the lattice constant I and the current velocity Cq = 1/t, the wave equation (H) becomes 



2cos6> 



2d 



,2 L 



+ \v2d+i(r)\ +\v2d+2(r) t / k=1 
2 (0(r,t + r) + 0(r,t-r)-20(r,t)) 

cos 8 



(Hr k ,t)~<j>(r,t)) 



2d+\v 2d+1 (r)\ - \v2d+i{r)\' 



2d+\v 2d +i{r)\ + |^d+2(r)r 



(*">*), 



where we have written [i = e l6 . 

Comparison between ([45]) and the Klein-Gordon equation (E 



provides 



c 2 (r) = c 2 



a 2 (r) = - 



2cos6» 



2d+ \v2d+i(r)\ + \ v2d+2(r) 



cos 9 



2d + I ^ 2 rf+i (r) I - I v 2 d+2 (r) | 
2rf+|^2d+i(r)| 2 + |^2d+2(r-)| 



(45) 



These results are an improvement over those obtained in section VII A 1 since a is a function of the position 
r. Moreover, c 2 (r) and a 2 (r) can vary independently from each other by properly adjusting the variations 
of V2d+ i(r) a nd V2d+2(r). Finally, the necessity of setting £2d+i = 1 which was found quite arbitrary in 



section VTIA, is no longer pertinent. The above results show that one really needs to keep both €2d+i 
and e2d+2 = — 1 simultaneously in order to obtain a complete version of the Klein-Gordon equation. 
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B. First order time equation (e e = — 1) 



The wave equation becomes 



^{(j)(r,t + T)-(j)(r,t-T)) = 
It 



I 2 I sin# \ ^ f 4>(r k , t) — (j>(r, t) 



t ^2d+|iA2d+i(r-)| 2 + |i/ 2 d +2 (r) 



E 



k=l 



2 ^ \ P 



sin6> / 2d + \v 2 d+i(r) \ 2 - |^ 2 d +2 (r) \ 2 
t V 2d+|^ 2(i+1 (r)| 2 + |j/ 2d+2 (r)| 2 , 

Comparison between ([l6|) and the Schrodinger equation (^0|) gives 

V(r) = _sinfl / 2d+\u 2 d+i{r)\ 2 -\v2d+2{r 
h T V2d+|^ 2rf+ i(r)| 2 + |^ d+2 (r)| 



r,t). (46) 




and 




sin# 



v 2 d+i{r)\ 2 + \ v 2 d+2{r)\ 2 



In contrast with the results discussed in section VII B, this new version of the discrete Schrodinger equation 



is now well defined. The potential V(r) and the mass m(r) can vary independently from each other by 
properly adjusting the variations of iy 2 d+i(f) and V2d+2(r). For instance, it is possible to formulate the 
Anderson model with diagonal and off-diagonal disorders separately. Of particular interest is the choice 
I ^2d+i(f ) | + | ^2d+2(f) | = constant which corresponds to keeping m(r) = constant while V{r) is varying 
with u 2d +i(r) and ^ 2 d +2 (r). 

Introducing two node currents could seem at first glance somehow arbitrary but is justified by the resulting 
well-defined Klein-Gordon and Schrodinger equations, equations (Ha) and (Ha) res pect ively. These equations 



do not suffer from the drawbacks of the 2d + 1 current model discussed in section VII 



IX. CONCLUSION 



It is at first glance a curiosity that starting from a discrete formulation of Huygens' principle and taking 
account of basic properties and symmetries, we have derived two equations that can be directly identified with 
the Klein-Gordon equation and the Schrodinger equation. On the one hand, this is not completely surprising 
considering that this simple model incorporates the basic features of wave propagation. On the other hand, the 
temporal and spatial symmetries which have been selected in building the model are simply those underlying 
the resulting discrete wave equations. The novelty of this approach is the systematic derivation resulting in 
a versatile unified equation, which properly tuned by a single parameter, yields two of the most fundamental 
wave equations in physics. So this derivation can be viewed as a general framework into which is merged the 
related approach of the TLM network model. In particular, it is worthwhile to point out that this unified 
equation can be derived on any discrete lattice or graph, thus underlining the generality of the construction 
(see Appendix). 

The most unusual feature of this formulation is the introduction of currents without any initial physical 
meaning. This can be compared with the cellular automata approach which is now routinely used for solving 
hydrodynamic problems. By starting from a minimal microscopic model of non physical particles which 
obey the most fundamental conservation laws of physics such as energy and momentum conservation, the 
macroscopic behavior of real fluids is recovered. The current model described in this paper has a similar 
status. A physical interpretation of the currents was given in the TLM network model by identifying them 
with electric field impulses traveling on a mesh of transmission lines 0]. Since in that particular case it 
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can be shown that the voltages and currents obey the same equations as the electric and magnetic fields 
of Maxwell's equations, it is not surprising to obtain Maxwell's equations from that model. However, such 
an interpretation fails in the abstract context we have considered in this paper, especially when we obtain 
the Schrodinger equation without any explicit analytic continuation from the diffusion equation (e.g. 1 i— > it 
p6|). Although they locally propagate as particles, and then, are scattered as waves, the currents cannot be 
interpreted either as fields or as particles. This is also in contrast with other approaches where the Schrodinger 
and Maxwell equations have been obtained from the dynamics of ensembles of Brownian particles 
thus providing a particle picture to standard wave equations. 

In considering further investigations, the natural question that comes immediatly to mind, concerns the 
possibility of describing other wave equations by the same approach. Some clues have already been indicated 
in this paper. For instance, we have discussed the possibility of describing higher order time derivatives of 
the field by releasing condition (|l3|). More generally, as the construction relies essentially on symmetries, 
we may wonder what equations would be obtained if a symmetry condition is not fulfilled or is replaced by 
another one. For instance, it will be easy to release the time-reversal symmetry and to solve numerically the 
Schrodinger equation within this assumption. This would be of fundamental interest in our understanding 
of strong localization in random media. Finally, there is no reason for limiting ourselves to scalar equations. 
As discussed in the introduction, the Maxwell equations as an example of vector wave equations, have been 
exhibited in the TLM approach. A further possibility is the description of spinor waves by introducing two 
different types of currents propagating on a Cartesian lattice. This work will be described in a forthcoming 
paper §. 
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APPENDIX: FORMULATION OF THE CURRENT-MODEL ON A GRAPH 

Although a n-dimensional Cartesian lattice has been considered throughout this paper, we have pointed 
out that most of our results did not depend on this particular lattice. This special kind of lattice, which 
involves an underlying manifold equipped with a metric, is required when considering the continuous limit of 
the discrete wave equation (Q). Indeed, in that limit, the lattice spacing I and the current velocity cq = 1/t 
are needed. In this appendix, we show that our model leads to the discrete wave equation (^4|) formulated on 
any arbitrary lattice. We are not going to reproduce in detail the derivation of Eq. (ff4|), which remains quite 
unchanged. Instead of this, we focus on the main differences between both approaches on either a Cartesian 
lattice or on any arbitrary lattice. 

We consider the most general lattice defined as a random lattice where each node is connected to z(r) 
neighbor nodes, z{r) denotes the coordination number. The notation z(r) indicates explicitly that the number 



of neighbors varies from node to node, as sketched in Fig. II. It is important to emphasize that no length 
is attributed to the bonds linking any two nodes. Thus, the bonds connected to one node represent only the 
connectivity of this arbitrary node. 
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FIG. 11. An arbitrary discrete lattice. 



As explained in the introduction, the currents propagate in one time step between a pair of nodes of a 
given bond. By taking into account the node currents, z(r) + 2 currents Ek(r,t) are incident on node r at 
each time step. Hence, a (z(r) + 2) x (z(r) + 2) scattering matrix is attached to node r. The definition (0) 
of the field becomes 



z(r)+2 

E 

k=l 



\ k {r)E k {r,t). 



(Al) 



The number of currents in the definition of the field depends now on the node r. Keeping in mind these 
definitions, and using the previous notation for the neighbors of node r, as described in Figure ^, it is eas; 
to check that all the steps of the derivation of equation (|4|) are left unchanged. For instance, conditions 
and ( |l3|) still hold to derive a closed field equation analogous to equation ([IB]). The new closed field equation 
reads straightforwardly 



z(r)+2 



k=l 



k=l 



ip(r,t-r). 



(A2) 



Then, the implementation of time-reversal symmetry and reciprocity can exactly be reproduced by substitut- 
ing each summation over 2d + 2 currents by a summation over z(r) + 2 currents. Isotropy is also introduced 
in the same way by noticing again that the number of pro pagating currents involved in this local symmetry 
now depends on the node r. As discussed in section [Vj, the lattice being arbitrary, equation (g9|) holds 
and the sign t\ can be fixed equal to 1. Then, we introduce naturally the ratios v z ^ +1 {r) = X z (r)+i/ x i 
and v z ( r ) +2 (r) = A z(r ) +2 /Ai which satisfy the conditions v* (r)+1 (r) / v z ( r ) + x(r) = e z {r)+i = ±1 and 



^V)+ 2 ( r )/^M+2( r ) 
context and read 



0(r,t) =Ai(r) 



±1. Finally, the equations (|42j), 



and 



are recovered in this 



z(r) 

E 

k=l 



Mr) 



E k (r, t) + v z{r ) +1 (r)E z{r}+1 (r, t) + v z ( r )+2{r)E z[r)+2 {r , t) 



z{r) + \v z(r)+ i{ r )\ 2 + K(r)+ 2 (r)r' 



r,t + r)+ e e 4>(r, t- r) = 

z(r) 

Ai(r) Y, W r ^t) - 4>(r,t)) + (z{r) + v 2 z{r)+l {r) + v 2 z[r)+2 (r)) <j>{r,t) 



k=l 



(A3) 



(A4) 



(A5) 
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where e e = ±1 and /i = e l9 . 

In equation ( A5), one recognizes the term J2k=l (4>{ r k, t) — <ft(r, t)) as a topological Laplacian defined on a 
graph. Thus, equation ( |A5j ) is a wave propagation equation formulated for a scalar field defined on a graph. 
This result can be compared to the attempts of formulating discretized field theories for scalar fields |2Q,M . 
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